前回課題の回答例

【問題】地震，プレート境界線，火山データのすべてを可視化してください

ヒント：

火山データの標高は可視化しなくて構いません
見やすくなるよう描画順や色等を工夫してください

回答例

課題に必要なコードだけで構成されていることを確認してください．

※関数の呼び出しがどのパッケージで行なわれているかを分かりやすくするため，パッケージ名::関数名という形でコーディングしてあります．通常はパッケージ名::の部分は省略して構いません（複数パッケージで関数名が重複する場合は必要です）．

#１．環境設定

#ggplot2パッケージはグラフ描画に必要です
if(!require(ggplot2)){
  install.packages("ggplot2")
  library(ggplot2)
}
# mapsパッケージは世界地図の描画に必要です
if(!require(maps)){
  install.packages("maps")
  library(maps)
}
# maptoolsパッケージはshapeファイルの読み込みに必要です
if(!require(maptools)){
  install.packages("maptools")
  library(maptools)
}
# broomパッケージはSpatialLinesDataFrameをdata.frameに変換するのに必要です
if(!require(broom)){
  install.packages("broom")
  library(broom)
}

#フォントファミリを設定する．（游ゴシック体）
if(.Platform$OS.type=="windows")
  # Windowsの場合
  windowsFonts(yugo=windowsFont("Yu Gothic"))
if(capabilities("aqua"))
  # MacOSの場合
  quartzFonts(yugo=quartzFont(rep("YuGo-Medium",4)))

#２．データの読み込み

#プレート境界線データの読み込み
# 発散型境界shapeファイルの読み込み 
ridge <-maptools::readShapeSpatial("plate/ridge.shp") 
# ggplot2での描画のためSpatialLinesDataFrameをdata.frameに変換する
ridge.fort<-broom::tidy(ridge)
# トランスフォーム型境界shapeファイルの読み込み 
transform <-maptools::readShapeSpatial("plate/transform.shp") 
# ggplot2での描画のためSpatialLinesDataFrameをdata.frameに変換する
transform.fort<-broom::tidy(transform)
# 収束型境界shapeファイルの読み込み
trench <-maptools::readShapeSpatial("plate/trench.shp") 
# ggplot2での描画のためSpatialLinesDataFrameをdata.frameに変換する
trench.fort<-broom::tidy(trench)

#火山データの読み込み
volc <- read.csv("GVP_Volcano_List.csv",header =TRUE,skip = 1,stringsAsFactors = FALSE)

#地震データの読み込み
quake<-read.csv("20110311JSTquery.csv",header = TRUE,stringsAsFactors = FALSE)

#３．描画

# ggplotオブジェクトを初期化する
g<-ggplot2::ggplot()

#プレート境界の描画
#発散型境界を描画する
g <- g + ggplot2::geom_path(data= ridge.fort, mapping = aes(x=long, y=lat, group=group),color="red")
#トランスフォーム型境界を描画する
g <- g + ggplot2::geom_path(data= transform.fort, aes(x=long, y=lat, group=group),color="red")
#収束型境界を描画する
g <- g + ggplot2::geom_path(data= trench.fort, aes(x=long, y=lat, group=group),color="red")

#火山の描画
g<-g+ggplot2::geom_point(data=volc,mapping = aes(x=Longitude,y=Latitude),shape=2,color="#000000")

#世界地図の描画
g<-g+ggplot2::borders(database = "world")

#地震の描画
g<-g+ggplot2::geom_point(data = quake,mapping = aes(x=longitude,y=latitude,color=depth,size=mag))

#使用するフォントファミリを設定する
g<-g+ggplot2::theme_bw(base_family = "yugo")

# 軸ラベルを描画する
g<-g+ggplot2::ylab("緯度")+ggplot2::xlab("経度")

#色の範囲を設定する．地震の発生深度が浅ければ赤，深ければ黄色とする．
g<-g+ggplot2::scale_color_continuous(low=rgb(1,0,0),high=rgb(1,1,0))

#点のサイズを0から6の範囲とする
g<-g+ggplot2::scale_size_continuous(range=c(0,6))

# タイトルを描画する
g <- g +ggplot2::ggtitle("プレート境界と地震と火山データ")

# ggplotオブジェクトを画面に描画する
g

#画像を保存する
ggsave("ex01.png",plot = g)

【発展問題】上記内容を特定の国に絞って可視化してください

ヒント：

ggplot2 の可視化データはすべてdata.frameです
subsetを使います

回答例

課題に必要なコードだけで構成されていることを確認してください．

#１．環境設定

#ggplot2パッケージはグラフ描画に必要です
if(!require(ggplot2)){
  install.packages("ggplot2")
  library(ggplot2)
}
# mapsパッケージは世界地図の描画に必要です
if(!require(maps)){
  install.packages("maps")
  library(maps)
}
# maptoolsパッケージはshapeファイルの読み込みに必要です
if(!require(maptools)){
  install.packages("maptools")
  library(maptools)
}
# broomパッケージはSpatialLinesDataFrameをdata.frameに変換するのに必要です
if(!require(broom)){
  install.packages("broom")
  library(broom)
}

#フォントファミリを設定する．（游ゴシック体）
if(.Platform$OS.type=="windows")
  # Windowsの場合
  windowsFonts(yugo=windowsFont("Yu Gothic"))
if(capabilities("aqua"))
  # MacOSの場合
  quartzFonts(yugo=quartzFont(rep("YuGo-Medium",4)))


#２．データの読み込み

#プレート境界線データの読み込み
# 発散型境界shapeファイルの読み込み 
ridge <-maptools::readShapeSpatial("plate/ridge.shp") 
# ggplot2での描画のためSpatialLinesDataFrameをdata.frameに変換する
ridge.fort<-broom::tidy(ridge)
# トランスフォーム型境界shapeファイルの読み込み 
transform <-maptools::readShapeSpatial("plate/transform.shp") 
# ggplot2での描画のためSpatialLinesDataFrameをdata.frameに変換する
transform.fort<-broom::tidy(transform)
# 収束型境界shapeファイルの読み込み
trench <-maptools::readShapeSpatial("plate/trench.shp") 
# ggplot2での描画のためSpatialLinesDataFrameをdata.frameに変換する
trench.fort<-broom::tidy(trench)

#火山データの読み込み
volc <- read.csv("GVP_Volcano_List.csv",header =TRUE,skip = 1,stringsAsFactors = FALSE)

#地震データの読み込み
quake<-read.csv("20110311JSTquery.csv",header = TRUE,stringsAsFactors = FALSE)


#３．データ加工

#日本地図を取得する
jp<-ggplot2::borders(database = "world",regions = "japan")
# 取得した日本地図の経度値の最大を取得する
max_x <- max(jp$data$long)
# 取得した日本地図の緯度値の最大を取得する
max_y <- max(jp$data$lat)
# 取得した日本地図の経度値の最小を取得する
min_x <- min(jp$data$long)
# 取得した日本地図の緯度値の最大を取得する
min_y <- min(jp$data$lat)

#地震データを日本地図の範囲内で発生したものに限定する
quake2<-subset(quake,(quake$longitude < max_x & quake$longitude > min_x)&(quake$latitude < max_y & quake$latitude > min_y))

#プレート境界線を日本地図の範囲内のものに限定する
ridge.fort2 <- subset(ridge.fort,(ridge.fort$long < max_x & ridge.fort$long > min_x)&(ridge.fort$lat < max_y & ridge.fort$lat > min_y))
transform.fort2 <- subset(transform.fort,(transform.fort$long < max_x & transform.fort$long > min_x)&(transform.fort$lat < max_y & transform.fort$lat > min_y))
trench.fort2 <- subset(trench.fort,(trench.fort$long < max_x & trench.fort$long > min_x)&(trench.fort$lat < max_y & trench.fort$lat > min_y))

#火山データを日本地図の範囲内のものに限定する
volc2 <- subset(volc,(volc$Longitude < max_x & volc$Longitude > min_x)&(volc$Latitude < max_y & volc$Latitude > min_y))

#４．描画

# ggplotオブジェクトを初期化する
g<-ggplot2::ggplot()

#プレート境界の描画
#発散型境界を描画する
g <- g + ggplot2::geom_path(data= ridge.fort2, mapping = aes(x=long, y=lat, group=group),color="red")
#トランスフォーム型境界を描画する
g <- g + ggplot2::geom_path(data= transform.fort2, aes(x=long, y=lat, group=group),color="red")
#収束型境界を描画する
g <- g + ggplot2::geom_path(data= trench.fort2, aes(x=long, y=lat, group=group),color="red")

#火山の描画
g<-g+ggplot2::geom_point(data=volc2,mapping = aes(x=Longitude,y=Latitude),shape=2,color="#000000")

#日本地図の描画
g<-g+jp

#地震の描画
g<-g+ggplot2::geom_point(data = quake2,mapping = aes(x=longitude,y=latitude,color=depth,size=mag))

#使用するフォントファミリを設定する
g<-g+ggplot2::theme_bw(base_family = "yugo")

# 軸ラベルを描画する
g<-g+ggplot2::ylab("緯度")+ggplot2::xlab("経度")

#色の範囲を設定する．地震の発生深度が浅ければ赤，深ければ黄色とする．
g<-g+ggplot2::scale_color_continuous(low=rgb(1,0,0),high=rgb(1,1,0))

#点のサイズを0から6の範囲とする
g<-g+ggplot2::scale_size_continuous(range=c(0,6))

# タイトルを描画する
g <- g +ggplot2::ggtitle("プレート境界と地震と火山データ")

# ggplotオブジェクトを画面に描画する
g

#画像を保存する
ggsave("ex02.png",plot = g)

地図情報の利用

前回は簡素な地図に観測データを描画しました．国単位での描画は可能でしたが，柔軟性には欠けていました．

本日は ggmap というライブラリを使って Google Map を用いた可視化を試みます．

はじめてのggmap

ggmap パッケージを用いると R で Google Map にプロットすることができるようになります．

それでは早速 Google Map を描画してみましょう．

#ライブラリを読み込む
if(!require(ggmap)){
  install.packages("ggmap")
  library(ggmap)
}

# 経度，緯度を指定する．
loc<- c(lon=139.4451841,lat=35.6944078)
# google map を取得する
map <- ggmap::get_map(location = loc)
# google map を描画する
ggmap::ggmap(map)

上記では（経度：139.4451841,緯度：35.6944078）を中心としたGoogleMapを描画しています．

下記のようなエラーが表示される場合は ggmap を再インストールすると修正されるかもしれません．

> ggmap(map)
 layer(mapping = NULL, data = NULL, stat = "identity", geom = <environment>,  でエラー: 
   使われていない引数 (geom_params = list(raster = c("#898C84", "#7D817D", "#818581", "#858981", "#8D9189", "#A6AAA2", "#B7BDB3", "#AFB3A7", "#B7BDB3", …

地図画像の取得

地図画像を取得するには get_map 関数を用います．引数を変更すれば地図の種類を変更することができます．

先程のコードを以下のように変更してください．変更箇所は get_map 関数の「maptype=“terrain”」とある場所だけです．

#ライブラリを読み込む
if(!require(ggmap)){
  install.packages("ggmap")
  library(ggmap)
}

# 経度，緯度を指定する．
loc<- c(lon=139.4451841,lat=35.6944078)
# ここを変更する．地形図を表示する
map <- ggmap::get_map(location = loc,maptype="terrain")
# google map を描画する
ggmap::ggmap(map)

同じようにmaptypeを変更するとGoogleMapの地図種類を変更できます．

航空写真としてプロットする

「maptype=“satellite”」で航空写真としてプロットします．

#航空写真 (Google Map)
map <- ggmap::get_map(location = loc,maptype="satellite")
ggmap::ggmap(map)

道路地図 (Google Map)としてプロットする

「maptype=“roadmap”」で道路地図としてプロットする．

#道路地図 (Google Map)
map <- ggmap::get_map(location = loc,maptype="roadmap")
ggmap::ggmap(map)

航空写真+道路地図 (Google Map)としてプロットする

「航空写真+道路地図」としてプロットします．

#航空写真+道路地図 (Google Map)
map <- ggmap::get_map(location = loc,maptype="hybrid")
ggmap::ggmap(map)

ズームレベルを変更する

「zoom=16」でズームレベルを「16」に変更します．

#ズームレベル16(ズームレベル1,２は使いものにならない)
map <- ggmap::get_map(location = loc,maptype="hybrid",zoom=16)
ggmap::ggmap(map)

モノクロにする

「color=“bw”」でモノクロとしてプロットします．

#モノクロ地図
map <- ggmap::get_map(location = loc,zoom=14,color="bw")
ggmap::ggmap(map)

OpenStreetMapから取得する．

「source = “osm”」でOpenStreetMap(http://www.openstreetmap.org/)から地図画像を取得してプロットします．

#モノクロ地図
map <- ggmap::get_map(location = loc,zoom=14,source = "osm")
ggmap::ggmap(map)

OpenStreetMap でも「zoom=16」や「color=“bw”」も指定することができます．

地図画像の取得2

マーカを設置する

以下の条件でマーカを設置した状態でGoogleMapを取得します．

説明	経度	緯度
中心点	139.4451841	35.6944078
マーカ1	139.4464	35.699
マーカ2	139.4454	35.695
マーカ3	139.4466	35.68142

loc<- c(lon=139.4451841,lat=35.6944078)
# マーカーを置く
df <- data.frame(lon=c(139.4464,139.4454,139.4466),lat=c(35.699,35.695,35.68142))
df

##        lon      lat
## 1 139.4464 35.69900
## 2 139.4454 35.69500
## 3 139.4466 35.68142

map <- ggmap::get_googlemap(center=loc,zoom=14,maptype="roadmap",markers=df)
ggmap::ggmap(map)

パスを書く

以下の条件でパスを描画した状態でGoogleMapを取得します．

説明	経度	緯度
中心点	139.4451841	35.6944078
点1	139.4464	35.699
点2	139.4454	35.695
点3	139.4466	35.68142

# パスを描く
df <- data.frame(x=c(139.4464,139.4454,139.4466),y=c(35.699,35.695,35.68142))
map <- ggmap::get_googlemap(center=loc,zoom=14,maptype="roadmap",path=df)
ggmap::ggmap(map)

緯度，経度を調べる

geocode 関数で位置（緯度，経度）を調べることができます．

loc<-ggmap::geocode(iconv("一橋大学",to="utf8"),output = "latlon",source = "google")
loc

##        lon     lat
## 1 139.4468 35.6948

geocode 関数は内部でGoogle Maps Geocoding API(https://developers.google.com/maps/documentation/geocoding)を呼び出しています．

手動で呼び出すには以下のようにします． URLを呼び出して戻り値としてのJSON形式のテキストを読み込み処理します．

if(!require(RCurl)){
  install.packages("RCurl")
  library(RCurl)
}
if(!require(jsonlite)){
  install.packages("jsonlite")
  library(jsonlite)
}
loc <- "一橋大学"
#文字コードを変換します
loc<-iconv(loc,to="UTF8")
#Google Maps Geocoding API URLを作ります
url<-paste("https://maps.googleapis.com/maps/api/geocode/json?address=",loc,sep = "")
url

## [1] "https://maps.googleapis.com/maps/api/geocode/json?address=一橋大学"

#URLエンコード
url<-URLencode(url)
url

## [1] "https://maps.googleapis.com/maps/api/geocode/json?address=%E4%B8%80%E6%A9%8B%E5%A4%A7%E5%AD%A6"

#APIを呼び出す
ret<-RCurl::getURL(url)
# JSONをパースする
ret<-jsonlite::fromJSON(ret)
# list形式をvector形式に変換後、(経度、緯度)となっているところを(緯度、経度)にrev関数を用いて変更します。
loc<-rev(unlist(ret$results$geometry$location))
loc

##      lng      lat 
## 139.4468  35.6948

get_openstreetmap 関数

OpenStreetMapに特化した get_openstreetmap 関数もありますが、Bounding Box（取得する地図の領域）指定が必要になるため少し不便です。

Bounding Boxは例えば以下のサイトで調べることができます。

http://boundingbox.klokantech.com/

map<-ggmap::get_openstreetmap(bbox=c(left=139.436728,bottom=35.689492,right=139.450986,top=35.700453),scale = OSM_scale_lookup(zoom = 16))
ggmap::ggmap(map)

ggplot2 と組み合わせる

ggmap は内部で ggplot2 を呼び出しているのでこれまでどおりのやり方でプロットできます．

loc<- c(lon=139.4451841,lat=35.6944078)
map <- ggmap::get_map(location =loc ,zoom=16,maptype="hybrid",source="google")
#変数にggmapの戻り値を代入する
g <- ggmap::ggmap(map)
#兼松講堂にプロットするためのデータフレーム
df <- data.frame(lon=c(139.44541814917307),lat=c(35.6949493737185),size=c(10))
#点をプロット。色は赤
g <- g + ggplot2::geom_point(data=df,aes(x=lon,y=lat,size=size,color="red"))
#今回は凡例の意味がないので消す
g + ggplot2::theme(legend.position = "none")

地図範囲を取得する

get_mapの戻り値の属性「bb」を参照すると，地図範囲を取得することができます．

loc<- c(lon=139.4451841,lat=35.6944078)
map <- ggmap::get_map(location =loc ,zoom=16,maptype="hybrid",source="google")
#地図範囲の取得
bb<-attr(map,"bb")
bb

##     ll.lat   ll.lon   ur.lat   ur.lon
## 1 35.68882 139.4383 35.69998 139.4521

この値を利用すれば地球規模の地震観測データから特定の領域のデータを取得することができるようになります．

地震観測データの可視化

#ライブラリを読み込む
if(!require(ggmap)){
  install.packages("ggmap")
  library(ggmap)
}

dat <- read.csv("20110311JSTquery.csv")
map <- ggmap::get_map(location=c(139.4452,35.69366),zoom=3,source="google")
bb<-attr(map,"bb")
dat2 <- subset(dat,latitude > bb$ll.lat & latitude < bb$ur.lat & longitude > bb$ll.lon & longitude < bb$ur.lon )
g <- ggmap::ggmap(map) 

g <- g + ggplot2::geom_point(data=dat2,aes(x=longitude,y=latitude,size=mag,color=depth),alpha=0.5)
g <- g + ggplot2::scale_color_continuous(low="#1F75FE",high="#00008B")
g <- g + ggplot2::scale_size_continuous(range=c(1,20))
g <- g + ggplot2::geom_text(data=dat2,aes(x=longitude,y=latitude,label=mag),size=2)
g

まとめ

本日は ggmap を用いた可視化手法を学びました．可視化結果自体は前回と変わりませんがGoogle Map を利用するだけで大部見映えがよくなったかと思います．

課題

　以下の問題を解いてください．その際に用いたRスクリプトファイルを提出してください．

※Rスクリプトファイルは可視化毎に作成してください．可視化に必要な内容のみを含めるようにしてください．

【問題】地震，プレート境界線，火山データのすべてをggmapを用いて可視化してください

ヒント：

地震観測データの可視化のコードを拡張することで実現できます
プレート境界線，火山データの subset は作る必要はありません．プロット時に自動的に削除されます．

回答例：

【発展問題】産総研の活断層データベースからデータを取得し，ggmap を利用して日本の活断層を描画してください．

ヒント：

データのダウンロード方法

活断層データベースにアクセスします．
「起震断層・活動セグメント検索」をクリックします．
検索方法で「すべてを表示」となっていることを確認し，「活動セグメント」をクリックします．
「検索結果をCSV形式でダウンロード」をクリックします．

read.csv読み込み後にデータを加工する必要があります．x座標は経度，y座標は緯度とする必要があります．
活断層の描画には geom_path を用います．geom_path の使い方は下記サンプルコードを参照してください．活断層データベースではサンプルコードのidと同じ役割として「活動セグメント番号」が使えます．

# ダミーデータを用意します
df<-data.frame(x=c(1,2,3,4),y=c(1,2,3,4),id=c(1,1,2,2))
df

##   x y id
## 1 1 1  1
## 2 2 2  1
## 3 3 3  2
## 4 4 4  2

# x,y だけ指定してgeom_pathで線を描画します
ggplot()+geom_path(data=df,aes(x=x,y=y))

# group を指定して描画します
ggplot()+geom_path(data=df,aes(x=x,y=y,group=id))

回答例：

戻る

Masaharu Hayashi を著作者とするこの作品はクリエイティブ・コモンズの表示 4.0 国際ライセンスで提供されています。

Day08: 地図情報を利用した可視化

前回課題の回答例

【問題】地震，プレート境界線，火山データのすべてを可視化してください

回答例

【発展問題】上記内容を特定の国に絞って可視化してください

回答例

地図情報の利用

はじめてのggmap

地図画像の取得

航空写真としてプロットする

道路地図 (Google Map)としてプロットする

航空写真+道路地図 (Google Map)としてプロットする

ズームレベルを変更する

モノクロにする

OpenStreetMapから取得する．

地図画像の取得2

マーカを設置する

パスを書く

緯度，経度を調べる

get_openstreetmap 関数

ggplot2 と組み合わせる

地図範囲を取得する

地震観測データの可視化

まとめ

課題

【問題】地震，プレート境界線，火山データのすべてをggmapを用いて可視化してください

【発展問題】産総研の活断層データベースからデータを取得し，ggmap を利用して日本の活断層を描画してください．